'''
Created on Jul 7, 2010

@author: oabalbin
'''
import numpy as np
from optparse import OptionParser
from datetime import datetime
from collections import defaultdict, deque

import RNAseq.arraytools.armanipulation as ca
import RNAseq.arraytools.subset_genes as sg
import RNAseq.array as seqarray

if __name__ == '__main__':
    
    optionparser = OptionParser("usage: %prog [options] ")
    optionparser.add_option("-f", "--annotFile", dest="annotFile",
                            help="annotation file for all files to use")
    optionparser.add_option("-d", "--outfolder", dest="outfolder",
                            help="output folder")
    optionparser.add_option("-s", "--normalsFile", dest="normalsFile",
                            help="annotation file with the normal samples")
    optionparser.add_option("-t", "--tissueName", dest="tissueName",
                            help="name of the tissue analyzed")

   
    (options, args) = optionparser.parse_args()
    
    
    t = datetime.now()
    tstamp = t.strftime("%Y_%m_%d_%H_%M")
    
    outfiles = ['outfile_expression_matrix']
    time_stamp_outfiles=[]
    for name in outfiles: 
        time_stamp_outfiles.append( options.outfolder+tstamp+'_'+options.tissueName+'_'+name+'.txt')
        
    thisfile = open(options.annotFile)
    outfile_matrix= open(time_stamp_outfiles[0],'w')
    
    #gene_list_folder = '/data/metadata/my_gene_lists/'
    #gene_list_folder = '/data/projects/prelims/genelists/'
    #gene_list_folder = '/data/projects/mirnas/noveltus/pol3/'
    gene_list_folder = '/data/projects/mirnas/noveltus/trnas'
    #gene_list_folder = '/data/projects/mirnas/noveltus/tuslist'
    
    #samples_list_file='/data/projects/prelims/lung_cell_lines.txt'
    #samples_list_file='/data/projects/prelims/lung_celllines_curated.txt'
    #samples_list_file='/data/projects/prelims/lung_cell_lines_rnaprot_matched.txt'
    #samples_list_file='/data/projects/prelims/lung_cell_lines_kras_localized.txt'
    # Phospho samples
    #samples_list_file='/data/projects/prelims/lung_cell_lines_rnaprot_matched_all.txt'
    #samples_list_file='/data/projects/prelims/lung_cell_lines_rnaprot_matched_rasonly.txt'
    #samples_list_file='/data/projects/prelims/lung_cell_lines_rnaprot_matched_rasonly_rnaseq.txt'
    samples_list_file='/data/projects/mirnas/noveltus/prostate_tissue_samples.dat'
    

    ####################
    ##### Parameters
    normal_samples_list = ca.list_of_names(open(options.normalsFile))
    
    expression_cutoff=0.0       # minimum median RPKM of expression
    base_expression_cutoff= 50  # In percentil. It is used to report the expression value of a particular gen in a sample. 
                                # The median usually can be used, however because in many cases median=1, It is maybe desireable to use for e.g the75%. 
    median_shift=1.0
    sdv_factor=0.25             #0.25
    not_report_isoforms=True   # Report gene outlier isoforms
    
    ########### Program
           
    # create the meme array
    filename=options.annotFile
    matrixfilename=options.annotFile+'memmap'
    header_cols=2
    matobj = seqarray.read_matrix(filename, matrixfilename, header_cols)
    all_samples, all_genes, all_gene_intervals = \
    seqarray.create_matobj_dicts(matobj.cols,matobj.sample_names)
    
    # Get all the gene files in the folder
    thosefiles = sg.read_files_folder(gene_list_folder,'.txt')
    # Get a subset of samples    
    thissampleslist = ca.list_of_names(open(samples_list_file))
    print thissampleslist
    sample_subset, sample_cols, sample_subset_rev = sg.get_sample_subset(thissampleslist, all_samples)
    #normal_sample_subset, normal_sample_cols, normal_sample_subset_rev = sg.get_sample_subset(normal_samples_list, all_samples)
    normal_sample_subset, normal_sample_cols, normal_sample_subset_rev = sg.get_sample_subset(normal_samples_list, sample_subset_rev)
        
    # 
    all_gene_list_names=deque()

    for gene_list_file in thosefiles:
        
        gene_list_name = gene_list_file.split('/')[-1]
        all_gene_list_names.append(gene_list_name)
        
        #print gene_list_file
        # subset of genes and samples to consider
        thisgenelist = ca.list_of_names(open(gene_list_file))
        print thisgenelist
        #print thisgenelist
        gene_subset, gene_rows, gene_subset_intervals, gene_subset_rows = \
         sg.get_gene_subset(thisgenelist, all_genes, all_gene_intervals) #list(set(all_genes.keys())) thisgenelist
        
        #gene_subset, gene_rows, gene_subset_intervals, gene_subset_rows = \
        #sg.get_gene_subset(list(set(all_genes.keys())), all_genes, all_gene_intervals) #list(set(all_genes.keys())) thisgenelist

        
        
        #sample_subset, sample_cols = sg.get_sample_subset(thissampleslist, all_samples)
        #normal_sample_subset, normal_sample_cols = sg.get_sample_subset(normal_samples_list, all_samples)
        
        # create a sub array  using the genelist and the samplelist
        
        expression_matrix = matobj.exprs[gene_rows, :]
        expression_matrix = expression_matrix[:,sample_cols]
        
        thisGeneArray = sg.create_thisGenearray(expression_matrix, gene_subset, sample_subset, \
                                                list(normal_sample_cols), gene_subset_intervals)
        # Possible thing to do to array
        #thisGeneArray.copa_norm_across_samples(median_shift)
        #thisGeneArray.copa_norm_in_sample()
        #thisGeneArray.fold_change_accross_samples()
        
        # Current used transformations
        thisGeneArray.log2_norm(median_shift)
        #thisGeneArray.center_scale_norm()
        thisGeneArray.normalize_exp_bynormals(normal_sample_cols)
        
        #thisGeneArray.nsaf_normalization()
        #thisGeneArray.fold_change_accross_samples_nsaf()
        #thisGeneArray.fold_change_accross_samples_mean()
        #thisGeneArray.log2_norm_nsaf(median_shift)
        ############# Print expression matrix after Copa normalization
        std_threshold=25
        #ca.print_expression_matrix(thisGeneArray, outfile_matrix)
        ca.print_expression_matrix_stdfiltered(thisGeneArray, outfile_matrix, std_threshold)
        #ca.print_prot_expression_matrix_stdfiltered(thisGeneArray, outfile_matrix, std_threshold)
        #ca.print_prot_expression_matrix_stdfiltered_nsaf(thisGeneArray, outfile_matrix, std_threshold)
        #ca.print_expression_matrix_isoforms_stdfiltered(thisGeneArray, gene_subset_rows, outfile_matrix, std_threshold)